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O ' Abstract 

o . 

CN . We develop perturbative expansions to obtain solutions for the initial-value problems of two im- 

Qh! portant reaction-diffusion systems, viz., the Fisher equation and the time-dependent Ginzburg-Landau 

(TDGL) equation. The starting point of our expansion is the corresponding singular-perturbation solu- 
tion. This approach transforms the solution of nonlinear reaction-diffusion equations into the solution of 
a hierarchy of linear equations. Our numerical results demonstrate that this hierarchy rapidly converges 
to the exact solution. 
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^ : 1 Introduction 

^ Many physical problems are described by nonlinear partial differential equations (pdes). In particular, 

much research interest has focused upon two classes of nonlinear pdes, viz., 
O . (i) soliton-bearing equations, which arise in the context of completely-integrable infinite-dimensional 
Hamiltonian systems [1,2]; and 

(ii) reaction-diffusion equations, which arise in the context of pattern-forming systems where local reac- 
^ tions are combined with spatial diffusion [3-6]. 

^ There is no general framework for obtaining the solution of the initial-value problem for an arbitrary 

O nonlinear pde. For the class of soliton-bearing equations in (i), powerful techniques like the inverse scat- 
^ ! tering transform and Backlund transformations enable the solution of these equations for arbitrary initial 
\ conditions [1,2]. Essentially, these methods reduce the problem of solution of a nonlinear soliton equa- 
^ tion to a sequence of linear equations. In a related context, a classic example of linearization is provided 
c3 by the Cole-Hopf transformation [7], which transforms the nonlinear Burgers' equation [8] into the linear 
diffusion equation. 

For the class of reaction-diffusion equations in (ii) above, there are as yet no systematic methods of 
linearization. These pdes have the general form: 

aiV^(f,t) = /(v^) + v2^, (1.1) 

where il){r^ t) is an order-parameter field, e.g., population density, chemical concentration, magnetization, 
which depends on space (r) and time (t). The order parameter may be either scalar or vector, depend- 
ing upon the number of variables which describe the physical system. The order parameter evolves in 
time due to a local reaction, described by the nonlinear term /(^), in conjunction with spatial diffusion. 
Reaction-diffusion equations are ubiquitous in pattern-forming systems, ranging from chemical and bio- 
logical physics to materials science and metallurgy [3-6]. 

An important example of Eq. (1.1) is the Fisher equation with /(^) = il) — -ip"^, which describes the 
growth and saturation of a species population [9,10]. Another important reaction-diffusion model is the 
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time-dependent Ginzburg-Landau (TDGL) equation with /{'if) = — [11,12], where 'if is an n- 

component vector, ^ = {i'i,i'2, •••! VaO- The case with n = 1 describes the phase ordering dynamics of a 
ferromagnet which has been suddenly quenched from the paramagnetic phase to the ferromagnetic phase. 
The TDGL equation with n = 2 describes phase ordering dynamics in superconductors, superfiuids and 
liquid crystals [12]. In this paper, we will investigate the perturbative linearization of the Fisher and TDGL 
equations for arbitrary initial conditions. 

This paper is organized as follows. Section 2 provides an overview of relevant analytical results, primar- 
ily in the context of the Fisher equation. Section 3 discusses our linearization scheme and presents detailed 
numerical results therefrom. Finally, Section 4 concludes this paper with a summary and discussion of our 
results. 



2 Overview of Analytical Results 

Let us first consider the Fisher equation in an infinite domain, which has the following form [9]: 

dttlj{f, t) = ^ - ^2 ^ V^^. (2.1) 

In general, Eq. (2.1) is supplemented with some arbitrary initial condition ^(r, 0). Typically, we are inter- 
ested in the case with ^ > 0, as the order parameter describes population density which cannot be negative. 
As a matter of fact, Eq. (2.1) is unstable for tjj <0. The homogeneous solution tjj* = Q is an unstable fixed 
point (FP) of the dynamics. Fluctuations about t/j* — diverge exponentially and saturate to the stable FP, 
ijj* = 1. 

There is no general solution available for the initial- value problem of Eq. (2.1). However, some im- 
portant analytical results are known for the case with dimensionality d = 1. Kolmogorov et al. [10] found 
that the Fisher equation has stable traveling- wave solutions 'ijj{x, t) =ip{x — vt) (called dines), which are 
domain walls with 

(a) velocity v >2, and ■?/'(— oo, t) — 1, '^{oo, t) — 0; or 

(b) velocity v < —2, and ■?/'(— cxo, t) = 0, ipioo, t) = 1. 

The qualitative forms of these solutions are easily obtained through a phase-portrait analysis, but the ex- 
plicit analytic forms are unknown. 

In more general work, Aronson and Weinberger [13] considered the d = 1 version of Eq. (1.1). They 
focused on functions f{ip) which satisfy the conditions 

/(V)>0 for ^e[0,l], 
/(O) = /(I) = 0, 

/'(O) > > /'(I). (2.2) 

These authors demonstrated that a broad class of initial conditions tlj{x, 0), with sharp interfaces, converge 
to a traveling- wave solution with a definite speed c*, which satisfies 



2[/'(0)]'/2 < ^* < 2^1/2^ 



L — sup 



V'e[o,i]. (2.3) 



A physical explanation of this velocity- selection principle has been formulated by various authors [14], and 
is often referred to as the marginal stability hypothesis. However, this approach has proven inadequate in 
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various applications. A more comprehensive understanding, based on structural-stability arguments, is due 

to Paquette et al. [15]. 

In the case of the Fisher equation, the Aronson- Weinberger result demonstrates that a large class of 
initial conditions converge to the cline solution with v — ±2, whose functional form is unknown as yet. An 
analytic form for a cline solution was first obtained by Ablowitz and Zeppetella [l6]forv — > 2. 

It is of obvious interest to obtain analytic forms for other cline solutions, particularly the case v = ±2. 
It is of even greater interest to obtain a general solution for the initial- value problem of Eq. (2.1). In 
this context, Puri et al. [17] have used singular-perturbation techniques, developed by Suzuki [18] and 
Kawasaki et al. (KYG) [19] in the context of the TDGL equation, to obtain an approximate solution for 
the initial-value problem of the Fisher equation: 

V'o(r,t) = 



V'L(f,t) = e*(^+^'V(r,0), (2.4) 

where 'tpL^f, t) is the solution of the linear part of the Fisher equation, and diverges with time. 

The approximate solution in Eq. (2.4) has a number of attractive features. For example, it is obtained for 
arbitrary dimensionality. Furthermore, initial conditions with sharp interfaces evolve into a traveling-wave 
front with asymptotic speed v — 2, in accordance with the Aronson- Weinberger result. The approach to 
the asymptotic velocity is as follows [17]: 

v{t) ^2-^^. (2.5) 

Unfortunately, this is not in agreement with the exact result for the d = 1 Fisher equation obtained by 
Bramson [20]: 

vit) ^ 2 - ^. (2.6) 

More generally, the fronts obtained using the singular-perturbation approximation are appreciably sharper 
than the exact result (obtained numerically) [17]. Furthermore, the solution in Eq. (2.4) is unable to ac- 
curately resolve front-front interactions [17]. An improved approximation has been proposed by Puri and 
Bray [21]. However, this is only valid for the less interesting class of initial conditions where the order pa- 
rameter ^ > everywhere, i.e., the evolving system is not characterized by the formation and interaction 
of fronts. 

Finally, it is also relevant to discuss the singular-perturbation solution for the TDGL equation: 

dtjjj{r, t) ^ ijj - jjj^ + V^V- (2.7) 

As mentioned earlier, the corresponding solution was obtained by KYG [19], who generalized a diagram- 
matic technique developed by Suzuki [18]. The approximate solution for the scalar TDGL equation is as 
follows: 



V^4f,t)=e*(i+^')V(r,0). (2.8) 



Puri and Roland [22] have obtained the singular-perturbation solution for the 2-component TDGL equation. 
This result has been generalized by Bray and Puri [23] and Puri [24] to obtain the time-dependent structure 
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factor for the n-component TDGL equation. As in the case of the Fisher equation, the approximate solution 
in Eq. (2.8) is characterized by fronts which are too sharp. Furthermore, the singular-perturbation solution 
is unable to properly resolve domain- wall interactions [17]. 

In this paper, we undertake a perturbative improvement of the singular-perturbation solution of the 
Fisher and TDGL equations. We develop a hierarchy of linear equations which rapidly converges to the 
exact solution. Our primary goal is methodological, viz., demonstrating the equivalence between nonlinear 
reaction-diffusion equations and a sequence of linear equations. 



3 Analytical and Numerical Results 
3.1 Fisher Equation 

The solution of the Fisher-equation with constant initial conditions suggests the following nonlinear trans- 
formation: 

l + 0(r,t) 

where > is an auxiliary field [12]. We confine ourselves to the physically interesting case with 1 > 
ijj >0. The corresponding pde satisfied by 0(f, t) is 

The singular-perturbation approximation is equivalent to dropping the nonlinear term on the right-hand-side 
(RHS) of Eq. (3.2). In that case, the solution (modified slightly from Eq. (2.4)) is 



1 - V'(r,0) 



(3.3) 



The approximate solution in Eq. (3.3) reduces to that in Eq. (2.4) in the limit of small ^(f, 0). However, in 
contrast to the earlier solution, Eq. (3.3) yields the exact solution in the homogeneous case 0) = V'(0)> 
viz., 

e*'0(O) 
= l-^(0)Ve-^(Q) - 

We will use the approximate solution in Eq. (3.3) as the starting point of a perturbative expansion, which 
yields a hierarchy of linear equations. First, notice that the exact pde obeyed by the solution V'o(^, t) is 

dtMr, = V'o - V'o + V'^o + (3.5) 

l-ipo 

We decompose the exact solution of the Fisher equation as ^ = x/jq + 9q, where 6*0 is the correction to the 
singular-perturbation solution. Then, the pde obeyed by 6'o(r, t) is 



dt9o{f, = (1 - 2V'o)^o - Ol + V% - (3.6) 
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We assume that 9o is small, and break it up as Oq — + 9i, where ipi solves the linear part of Eq. (3.6): 

(1 + 00)-^ 

The exact pde obeyed by 9i (r , t) is then 

dMr, t)^[l- 2(V'o + i^iWi -9l + V% - i^l (3.8) 

which has the same general form as Eq. (3.6). 

This process can be continued indefinitely. At the n*^ level, we decompose as 6'„_i = 4>n + 9n, where 
i/jn solves the linear version of the pde obeyed by ^n-i- The corresponding pde obeyed by 9n is then 



n 



dt9n{f, t)=il-2j2^k]0n-9l + V'9n - (3.9) 

V fc=0 / 

The approximate perturbative solution at the n}^ level of iteration is then obtained as 

^(r,t) ^^^fc(r,t). (3.10) 

fc=0 

The general form of the pdes obeyed by V'fe (for k > 1) is 

dtMr, t) = a{f, t)il^k + VVfe + Kr, t), (3.11) 

where a(r, t) and 6(r, t) are functions of space and time. A formal solution of Eq. (3.11) reads (somewhat 
symbolically): 

V'fe = (^t - - a)~' & . (3.12) 

It can be expanded in a Taylor-series in a. Denoting R{x, t) as the functional inverse of {dt — V^), such 
that {dt - V^)i?(f , t) = 5{x)6{t), Eq. (3.12) can be written as 





t 



.Jm'I ,f' I ,f I ,r-Rif- fHr", m(r" - f", t' - fW, O + • • • 



(3.13) 

Note that for 6(f, t) = 5{t) 4>{f, 0) and a = 1, we recover the solution to the linear part of Eq. (3.2). 

We will demonstrate shortly that this perturbative expansion converges very rapidly to the exact solu- 
tion. However, we stress that our interest in such an expansion is more methodological than operational, 
i.e., the above procedure converts the problem of solution of the nonlinear Fisher equation to a hierarchy 
of linear equations. Of course, the same procedure would serve to improve any approximate solution - the 
singular-perturbation result is merely a convenient starting point. 
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3.2 Time-Dependent Ginzburg- Landau Equation 



Before we present our numerical results, it is relevant to discuss the corresponding procedure for the scalar 
TDGL equation. (The generalization to the n-component case is relatively straightforward.) The corre- 
sponding nonlinear transformation is obtained by solving the TDGL-equation for constant initial conditions 
and reads 



o 



(r. t) 



(3.14) 



where is the appropriate auxiliary field. In this case, we confine ourselves to the physically interesting 
case with |^| < 1. The pde satisfied by 0(r, t) is [12] 



2, 30(V0)2 



1 + 



(3.15) 



As before, the singular-perturbation solution is obtained by neglecting the nonlinear term on the RHS of 
Eq. (3.15) as follows: 



Vl + 0o(r,t)2' 



0o(f,t) = e*(i+^^) 



Vl-V'(r,0)2 



(3.16) 



This solution constitutes an improvement over the solution in Eq. (2.8) as it is exact in the homogeneous 
case. The exact pde satisfied by ipo is then 



3^o(V^o) 



We decompose the solution of the TDGL equation as V' = V'o + do, where the pde obeyed by 9o{f, t) is 

3^o(VV^o)' 



(3.17) 



(3.18) 



As before, we designate the solution of the linear part of this equation as "01 t)- If we decompose the 
overall solution as Oq = ijji + 6i, the pde obeyed by 9i (r, t) is 



dt9i{r, t)^[l- 3(V'o + - 3(V'o + tpiWl -9l + V% - (SV'o + V'OV'?- 



(3.19) 



Again, we can develop an infinite hierarchy of equations. The pde for 9n{r,t) at the n*^ level of this 
hierarchy is 



dt9n{r,t) 



9n-3[J2^k]9l-9^^ + V^9r, - 3^^^, + (3_20) 



.A;=0 



fc=0 



Therefore, the problem of solution of the nonlinear TDGL equation is again reducible to the solution of an 
infinite hierarchy of linear equations for ipk{^, t). The approximate perturbative solution at the n*^ level is 
given by Eq. (3.10). 
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3.3 Numerical Results for the Fisher Equation 



In this subsection, we will numerically examine the convergence properties of the hierarchy of linear equa- 
tions presented above. For the sake of brevity, we confine ourselves to presenting numerical results for the 
Fisher equation in d = 1, 2. Similar results are obtained for the TDGL equation also. 

Our numerical results for the Fisher equation, referred to as the "exact solution" subsequently, were 
obtained by implementing an Euler-discretized version of Eq. (2.1) (with an isotropic Laplacian) in d — 
1,2. In both cases, we used periodic boundary conditions. The discretization mesh sizes in d — 1 were 
At = 0.001 and Ax = 0.1 in time and space, respectively. The corresponding mesh sizes in d = 2 were 
At = 0.001 and Ax = 0.2. The lattice size was A^i = 40000 in d = 1, and = 1000^ in = 2. The 
spatial coordinates are x G [—2000,2000] in d = 1; and x,y E [—100, 100] in d = 2. Our perturbative 
solutions for 1^0,1^1,1^2, etc. were also obtained numerically by solving the relevant linearized equations. 
The discretization meshes and lattice sizes for the perturbative solutions are identical to those described 
above. 

Figure 1(a) shows the profile of a front arising from a seed initial condition {■'tp{x, 0) = 0.055(x)) for the 
d = 1 Fisher equation. The evolution gives two equivalent clines moving in opposite directions - we focus 
on the cline with t> > 0. In Figure 1(a), we show the exact solution (solid line); the singular-perturbation or 
n — solution (dashed line); the n = 1 solution (dotted line); and the n = 2 solution (dot-dashed line). The 
n — 3 perturbative solution (dot-dot-dashed line) is already numerically indistinguishable from the exact 
solution on the scale of the figure - we will quantify the error shortly. We have numerically confirmed that 
the perturbation series in Eq. (3.10) is strongly convergent, and the inclusion of higher-order terms (n > 3) 
does not change the solution appreciably. 

Let us next examine the time-dependence of the front velocity. Recall that the n = solution did not 
exhibit the correct approach to the asymptotic velocity v = 2. Figure 1(b) plots v{t) vs. for the solutions 
depicted in Figure 1(a). Again, we see that the n = 3 result is almost coincident with the exact result. 

Figures 2(a)-(b) study the evolution of the d = 1 Fisher equation from a random initial condition, 
consisting of 10 randomly-distributed seeds of random height (between and 0.1). In Figures 2(a)-(b), we 
focus on the collision and merger of two fronts, showing the exact solution, and the solutions for n = 0, 1 , 3. 
As in Figure 1(a), the n = 3 solution is numerically indistinguishable from the exact solution. 

To quantify the error involved in our approximations, we compute the "distance" between the exact 
solution ijjeix, t) and an approximate solution ilJa{x, t) as follows: 



where L is the lattice length. Figure 2(c) plots D(t) vs. t on a semi-logarithmic scale for the different 
solutions depicted in Figures 2(a)-(b). We obtain D{t) by averaging over 25 independent initial conditions 
constructed as the random superposition of seeds, Gaussians, sine-cosine functions, etc. The maximum 
error for the n = 3 solution is three orders of magnitude smaller than that for the n = (singular- 
perturbation) solution. This quantifies the rapid convergence of our perturbative hierarchy. Of course, 
the error asymptotically approaches zero for all solutions as ijj ^ ijj* = 1 everywhere. 

Next, we consider results for the d = 2 Fisher equation. Figure 3(a) shows an evolution snapshot 
obtained from a seed initial condition ip{x, y, 0) = 0.05S{x)S{y). The dark circular region refers to the 
exact solution, and denotes points where ip > 0.5. The solid line refers to the front position for the n = 
solution, and is defined by points where ip = 0.5. Figure 3(b) shows the corresponding variation of the order 
parameter along a horizontal cross-section of the lattice. As in the d = 1 case, the singular-perturbation 



L/2 




(3.21) 



-L/2 
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front is appreciably sharper than the exact front. For clarity, we do not show the higher-order perturbation 
results - as before, the n = 3 result is indistinguishable from the exact solution. We focus on the front 
position so as to clarify the convergence of the perturbative solution. Figure 3(c) plots the front velocity 
v{t) vs. for the solutions depicted in Figure 3(b). 

Finally, Figures 4(a)-(b) study the evolution from a random initial condition for the d — 2 Fisher 
equation. The initial condition consists of 25 randomly-distributed seeds with amplitudes between and 
0.1. Figure 4(a) shows an evolution snapshot at t = 12. (For clarity, we only show a 500^ corner of 
the 1000^ lattice.) The coding is the same as that for Figure 3(a). Figure 4(b) shows the variation of the 
order parameter along a horizontal cross-section of the snapshot shown partially in Figure 4(a). We present 
results for the exact solution, and approximate solutions with n — 0, 1, 3. The approximation error is 
quantified in Figure 4(c), where we plot D{t) vs. t on a semi-logarithmic scale. The error is obtained as 
the d = 2 generalization of the quantity defined in Eq. (3.21). As in the d = 1 case, the error was obtained 
as an average over 25 independent initial conditions of different types. 

4 Summary and Discussion 

Let us conclude this paper with a summary and discussion of the results presented here. We have studied 
two important examples of reaction-diffusion systems, viz., the Fisher equation and the time-dependent 
Ginzburg-Landau (TDGL) equation. We are interested in the linearization of these and other reaction- 
diffusion equations, i.e., conversion of the nonlinear problem to a linear problem. In both cases, we find 
that the singular-perturbation solution is a good starting point for a perturbative expansion. This expan- 
sion transforms the problem of solution of the nonlinear partial differential equation to the solution of 
a hierarchy of linear partial differential equations. Our numerical studies demonstrate that this hierarchy 
rapidly converges to the exact solution of the relevant equation. However, we should stress that our primary 
interests in this paper are methodological rather than operational. 

The techniques developed here are of general applicability to a wide range of reaction-diffusion equa- 
tions. Of course, any approximate solution to the initial- value problem for a given equation is a good 
starting point for a perturbative expansion. However, we find that the singular-perturbation solution ap- 
pears to be particularly convenient in that it yields an extremely accurate solution within a few steps of the 
perturbation expansion. 
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Figure 1: (a) Evolution snapshot of the d = 1 Fisher equation for a seed initial condition ip{x, 0) = 0.05(5(x). We plot the order 
parameter tp{x,t = 200) vs. x for the front moving with velocity v > 0. Details of our simulation are provided in the text. We 
present results for the exact solution (solid line); the singular-perturbation orn = solution (dashed line); the n — 1 solution 
(dotted line); the n = 2 solution (dot-dashed line); and the n = 3 solution (dot-dot-dashed line), (b) Plot of front velocity v{t) 
vs. for the solutions depicted in (a). The line-type usage is the same as that in (a). 
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Figure 2: (a) Evolution snapshot of the = 1 Fisher equation for a random initial condition, consisting of 10 randomly- 
distributed seeds of amplitude between and 0. 1 . We plot the order parameter iJj{x, t) vs. x att = 40, focusing on a front-front 
collision. The results shown are analogous to those in Figure 1(a), except we do not show the case n = 2. (b) Analogous to (a), 
but at the later time t = 80. (c) Plot of D{t) vs. f on a semi-logarithmic scale, where D{t) is the "distance" (defined in Eq. (3.21)) 
between the exact solution and an approximate solution. The error D{t) is obtained as an average over 25 independent initial 
conditions of various types. We show results for the approximate solutions with n = 0, 1, 3, using the same hne-types as in 
(a)-(b). 
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Figure 3: (a) Evolution snapshot of the d = 2 Fisher equation from a seed initial condition V'(.T, y,0) = 0.05S{x)S{y). The 
dark region corresponds to the exact solution alt = 40, and denotes points where tp > 0.5. The soUd Une refers to the singular- 
perturbation (n = 0) solution att = 40, and denotes points where V = 0.5. (b) Variation of the order parameter along a 
horizontal cross-section for the evolution depicted in (a). We plot t/j{x,y,t = 40) vs. x for y = 0, and focus on the front 
solution. The results presented are analogous to those in Figure 2(a). (c) Plot of front velocity v{t) vs. for the solutions 
depicted in (b). 
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Figure 4: (a) Evolution snapshot of the d = 2 Fisher equation from a random initial condition, consisting of 25 randomly- 
distributed seeds with amplitude between and 0.1. The dark regions denote the exact solution at t = 12, and the solid 
line denotes the corresponding singular-perturbation solution. For clarity, we only show a 500^ corner of the 1000^ lattice, (b) 
Variation of the order parameter along a horizontal cross-section of the snapshot partly shown in (a). We plot tp{x,y — 0,t — 12) 
vs. X for the entire range of x-values. The line-type usage is the same as earlier, (c) Analogous to Figure 2(c), but for the d = 2 
Fisher equation. 
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